#pragma once

using namespace System;
using namespace System::IO;
using namespace NUnit::Framework;
using namespace LatoolNet;

namespace LatoolNetTest {

	[TestFixture]
	public ref class EigenvalueProblemTest {
	public:

		[Test]
		void TestSymmetricEigenProblem() {

			int rownum = 4;
			int colnum = 4;
			Matrix^ a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);

			a[0, 0] = 1;
			a[0, 1] = 2;
			a[0, 2] = 3;
			a[0, 3] = 4;

			a[1, 1] = 2;
			a[1, 2] = 3;
			a[1, 3] = 4;

			a[2, 2] = 3;
			a[2, 3] = 4;

			a[3, 3] = 4;

			Matrix ^ orig = a->Clone();

			EigenvalueProblem ^ evp = gcnew EigenvalueProblem();
			evp->Solve(a);

			int numberOfLambda = evp->EigenValue->RowNum;

			for (int n = 0; n < numberOfLambda; n++) {
				double lambda = evp->EigenValue[n, 0];
				Matrix^ vector = gcnew Matrix(4, 1);
				for (int i = 0; i < 4; i++) {
					vector[i, 0] = evp->EigenVector[i, n];
				}

				Matrix^ lhs = orig * vector;
				Matrix^ rhs = lambda * vector;

				for (int i = 0; i < 4; i++) {
					Assert::AreEqual(lhs[i, 0], rhs[i, 0], 1e-10, "Test: EigenvalueProblem");
				}
			}
		};

		[Test]
		void TestGenralizedSymmetricDefiniteEigenProblem() {

			int rownum = 4;
			int colnum = 4;
			Matrix^ a = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);

			a[0, 0] = 0.24;
			a[0, 1] = 0.39;
			a[0, 2] = 0.42;
			a[0, 3] = -0.16;

			a[1, 1] = -0.11;
			a[1, 2] = 0.79;
			a[1, 3] = 0.63;

			a[2, 2] = -0.25;
			a[2, 3] = 0.48;

			a[3, 3] = -0.03;

			Matrix^ b = gcnew Matrix(rownum, colnum, MatrixType::DoubleSymmetric);
			b[0, 0] = 4.16;
			b[0, 1] = -3.12;
			b[0, 2] = 0.56;
			b[0, 3] = -0.10;

			b[1, 1] = 5.03;
			b[1, 2] = -0.83;
			b[1, 3] = 1.09;

			b[2, 2] = 0.76;
			b[2, 3] = 0.34;

			b[3, 3] = 1.18;

			Matrix ^ origa = a->Clone();
			Matrix ^ origb = b->Clone();

			EigenvalueProblem ^ evp = gcnew EigenvalueProblem();
			evp->GeneralizedSolve(a, b);

			int numberOfLambda = evp->EigenValue->RowNum;

			for (int n = 0; n < numberOfLambda; n++) {
				double lambda = evp->EigenValue[n, 0];
				Matrix^ vector = gcnew Matrix(4, 1);
				for (int i = 0; i < 4; i++) {
					vector[i, 0] = evp->EigenVector[i, n];
				}

				Matrix^ lhs = origa * vector;
				Matrix^ rhs = lambda * origb * vector;

				for (int i = 0; i < 4; i++) {
					Assert::AreEqual(lhs[i, 0], rhs[i, 0], 1e-10, "Test: EigenvalueProblem");
				}
			}
		};

	};
}
